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Abstract 

An 8 -node decohesion element implementing different 
criteria to predict delamination growth under mixed-mode 
loading is proposed. The element is used at the interface 
between solid finite elements to model the initiation and 
propagation of delamination. A single displacement-based 
damage parameter is used in a softening law to track the 
damage state of the interface. The power law criterion and a 
three-parameter mixed-mode criterion are used to predict 
delamination growth. The accuracy of the predictions is 
evaluated in single mode delamination and in the mixed- 
mode bending tests. 

Introduction 

Delamination is one of the predominant forms of 
failure in laminated composites due to the lack of 
reinforcement in the thickness direction. Delamination as a 
result of impact or a manufacturing defect can cause a 
significant reduction in the compressive load-carrying 
capacity of a structure. The stress gradients that occur near 
geometric discontinuities such as ply drop-offs, stiffener 
terminations and flanges (Figure 1), bonded and bolted 
joints, and access holes promote delamination initiation, 
trigger intraply damage mechanisms, and cause a 
significant loss of structural integrity. 

The fracture process of high performance composite 
laminates is quite complex, involving not only interlaminar 
damage (delamination), but also intralaminar damage 
mechanisms (e.g. matrix cracking, fiber fracture). 
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Figure 1. Experiment illustrating stiffener-flange 
debonding. 

The objective of this paper is to present a method to 
simulate progressive delamination with decohesion 
elements based on a new mixed-mode failure criterion. The 
numerical predictions are compared with analytical closed- 
form solutions and experimental results. 

Numerical Simulation of Delamination 

The study of delamination mechanics may be divided 
into the study of delamination initiation and the analysis of 
delamination propagation. Delamination initiation analyses 
are usually based on stresses and use criteria such as the 
quadratic interaction of the interlaminar stresses in 
conjunction with a characteristic distance 1, 2 . The 
characteristic distance is an averaging length that is a 
function of geometry and material properties, so its 
determination always requires extensive testing. 

Most analyses of delamination growth apply a fracture 
mechanics approach and evaluate strain energy release rates 
G for self-similar delamination growth. The G values are 
usually evaluated using the virtual crack closure technique 
(VCCT) proposed by Rybicki and Kanninen’. The VCCT 
technique is based on Irwin’s assumption that when a crack 
extends by a small amount, the energy released in the 
process is equal to the work required to close the crack to 
its original length. The Mode I, Mode II, and Mode III 
energy release rates can then be computed from the nodal 
forces and displacements obtained from the solution of a 
finite element model. The approach is computationally 
effective since the energy release rates can be obtained 
from only one analysis. 


In the present paper, an approach is proposed that is 
well suited to nonlinear progressive failure analyses where 
both ply damage and delaminations are present. The 
approach consists of placing interfacial decohesion 
elements between composite layers. A decohesion failure 
criterion that combines aspects of strength-based analysis 
and Fracture Mechanics is used to simulate debonding by 
softening the element. The proposed constitutive equations 
for the interface are phenomenological mechanical relations 
between the tractions and interfacial separations. With 
increasing interfacial separation, the tractions across the 
interface reach a maximum, decrease, and vanish when 
complete decohesion occurs. The work of normal and 
tangential separation can be related to the critical values of 
energy release rates 4 . 

Decohesion can be implemented as a material response 
such as the Dudgale-Barenblatt (D-B) type cohesive zone 5 . 
The D-B cohesive zone model was first applied to the 
analysis of concrete cracking by Hillerborg et al 6 . The 
concept has also been used by Needleman 7 to simulate fast 
crack growth in brittle solids. Needleman considered that 
cohesive zone models are particularly attractive when 
interfacial strengths are relatively weak when compared 
with the adjoining material, as is the case in composite 
laminates. 

In order to predict the initiation and growth of 
delaminations, an 8-node decohesion element shown in Fig. 
2 was developed and implemented in the ABAQUS finite 
element code. The development of this element is based on 
prior work 8, l) . The decohesion element is used to model the 
interface between sublaminates or between two bonded 
components. The element consists of a zero-thickness 
volumetric element in which the inteipolating shape 
functions for the top and bottom faces are compatible with 
the kinematics of the elements that are being connected to 
it. The material response built into the element represents 
damage using a cohesive zone ahead of crack tip to predict 
delamination growth. The concept of interface elements has 
been used in different types of problems: compression after 
impact 10, n , damage growth from discontinuous plies 12 , and 
diametrical compression of composite cylinders 1 ’. 



Figure 2. Eight-node decohesion element; t=0. 


Constitutive Equations 

The need for an appropriate constitutive equation in 
the formulation of the interface element is fundamental for 


an accurate simulation of the interlaminar cracking process. 
A constitutive equation is used to relate the traction cr to the 
relative displacement 8 at the interface. Some softening 
models that have been proposed are shown in Fig. 3 and 
include: linear elastic-perfectly plastic; linear elastic-linear 
softening; linear elastic -progressive softening; linear 
elastic-regressive softening; and Needleman 7 . One 
characteristic of all softening models is that the cohesive 
zone can still transfer load after the onset of damage ( 8 ° in 
Figure 3). For pure Mode I, II or III loading, after the 
interfacial normal or shear stresses attain their respective 
interlaminar tensile or shear strengths, the stiffnesses are 
gradually reduced to zero. The area under the stress-relative 
displacement curves is the respective (Mode I, II or III) 
fracture energy. Using the definition of the J integral 
proposed by Rice 14 , it can be shown that for small cohesive 
zones, 

a{S)dS = G c (1) 

where G c is the critical energy release rate for a particular 
mode, and S F is the corresponding relative displacement at 
failure (S pp , 8 pro , S, in S Ne , or 8 re in Figure 3). 



Figure 3. Strain softening constitutive models. 


Bilinear Softenins Model 

The linear elastic-linear softening (bilinear) model is 
the simplest to implement, and is most commonly used 10, 1 ’ 
151 . The double cantilever beam test shown in Figure 4 
illustrates the material response. Point 1 in Figure 4 is 
subjected to a low tensile load that is within the linear 
elastic range. A high initial stiffness K p (penalty) holds the 
top and bottom faces of the interface element together. 
Point 2 represents the onset of damage. In single -mode 
delamination, the traction at point 2 is equal to the 
corresponding interlaminar strength of the material, O c . As 
the relative displacement increases, the interface 
accumulates damage and the traction is lower than the 
strength (point 3). The energy released at point 3 is the area 
of the triangle 0-2-3. If the load were to reverse, point 3 
would unload to the origin, as shown in the figure. 

The critical value of the energy release rate is attained 
at point 4. For any relative displacement larger than point 4, 



(2) 


the interface does not carry any tensile or shear loads (point 
5). In other words, at point 4 all the available interfacial 
fracture energy has been completely consumed. Note that 
when modeling delamination with a softening response, the 
delamination tip is not defined explicitly. While the onset 
of damage occurs at point 2 in Figure 4, the delamination 
tip could be defined as the point where the tractions at the 
interface are zero, which is point 4. 

The softening response illustrated in Figure 4 is 
representative of the tension or the shear response but not 
compression. It is assumed that compression loads do not 
cause delamination or softening, and the effect of 
compression on damage of the interface was neglected in 
the present work. 



The problem of contact of the crack faces after failure 
is addressed by re-applying the normal penalty stiffness. 
The process of reapplying the normal stiffness when 
interpenetration is detected is typical of solution procedures 
of contact problems using penalty methods in a constrained 
variational formulation. 

The concept of decohesion zones to simulate 
delamination growth in composites is usually implemented 
by means of interface elements connecting the individual 
plies of a composite laminate. Decohesion elements can 
model the discontinuity introduced by the growth of 
delaminations. They can be divided into two main groups: 
continuous interface elements and point interface elements. 
Several types of continuous interface elements have been 
proposed, ranging from plane interface elements with zero 
thickness connecting solid elements 10 ’ 1 ; plane interface 
elements with finite thickness connecting shell elements 1 ’; 
line interface elements 12, 17 ’ 18 ; and spring interface elements 
that connect pairs of nodes 15, 16 . 

The bilinear interfacial constitutive response shown in 
Figure 4 can be implemented as follows: 

i) 8 < 8° => the constitutive equation is given by: 


a = K P S 

ii) 8° <8< 8 F => the constitutive equation is given by: 

a = {\-D)K p S (3) 

where D represents the damage accumulated at the 
interface, which is zero initially, and reaches 1 when 
the material is fully damaged. 

iii) 8 > 8 F => all the penalty stiffnesses are set equal to 
zero. If crack closure is detected, interpenetration is 
prevented by reapplying only the normal stiffness. 
Frictional effects are neglected. 

The properties required to define the bilinear 
interfacial softening behavior are the initial stiffness 
(penalty) K P , the fracture energies G IC , G I1C , and G mc and 
the corresponding nominal interlaminar tensile and shear 
strengths, T and S. The accuracy of the analysis depends on 
the penalty stiffness K P that is chosen. High values of K P 
avoid interpenetration of the crack faces but can lead to 
numerical problems. Several values have been proposed for 
the penalty stiffness, K P : 10 7 N/mm 3 [Ref. 10], 5.7x 10 7 
N/mm’ [Ref. 19], 10 8 N/mm 3 [Ref. 16], Other authors have 
determined the value for the penalty stiffness as a function 
of the interface properties. Daudeville et al. 2(> have modeled 
the interface as a resin rich zone of small thickness, t h and 
proposed a penalty stiffnesses defined as: 


tsI -E 3 . v u _ 2Gf 3 . v m _ 2G 2i 
K P —— ,A P — - ,Kp — - 


(4) 


where G 23 , G I3 , and E 3 are the elastic moduli of the resin 
rich zone. After a substantia] number of numerical 
experiments, Moura 21 detennined that a penalty stiffnesses 
of only 1 0 6 N/mm 3 for all modes produces essentially the 
same results while avoiding potential convergence 
problems during the nonlinear procedure. 

In order to fully define the interfacial behavior, the 
unloading response must be specified. Petrossian et al. 16 
have proposed an unloading curve with a slope 
corresponding to Hooke's law. Such a procedure, typically 
used in the formulation of plasticity problems, would lead 
to the use of the same penalty stiffness when reloading and 
to permanent relative displacements along the interface 
when the load reverts to zero. Crisfield et al. 17, 18 and 
Daudeville 20 , on the other hand, have proposed that with 
reversing strains the material unloads directly toward the 
origin, as shown in Figure 4. The assumption is that during 
reloading the interfacial stiffness is lower than the original 
(undamaged) stiffness. Such a procedure simulates the 
effects of the previous damage mechanisms that occurred 
along the interface and was therefore adopted in the present 
work. 



Mixed Mode Delamination Criterion 


(10) 


In structural applications of composites, delamination 
growth is likely to occur under mixed-mode loading. 
Therefore, a general formulation for interface elements 
must deal with mixed-mode delamination growth problems. 

Under pure Mode I, II or III loading, the onset of 
damage at the interface can be determined simply by 
comparing the stress components with their respective 
allowables. Note that the onset of damage does not imply 
the initiation of delamination, since the tractions closing the 
crack at onset are at their maximum value. Under mixed- 
mode loading, however, damage onset may occur before 
any of the stress components involved reach their 
respective allowables. 

A mixed-mode criterion is proposed here that is based 
on a few simplifying assumptions. Firstly, it is assumed that 
delamination initiation can be predicted using the quadratic 
failure criterion 




for a. < 0 


(5) 

(6) 


where O- is the transverse normal traction and T x: and T y:: are 
the transverse tractions. T and S are the nominal normal 
tensile and shear strengths, respectively. 

The delamination mechanisms in Mode II and Mode 
III are assumed to be the same. Therefore, Mode III can be 
combined with Mode II by using a total tangential 
displacement Su defined as the norm of the two orthogonal 
tangential relative displacements 8 and Sy as 

$n =pl+$l (7) 

The total mixed-mode relative displacement 8 m is defined as 

= V <5- 2 + 8u (8) 

where 8 : is the relative opening (Mode I) displacement. 
Using the same penalty stiffness in Modes I and II, the 
tractions are 

<? : =K p 8 : 

t xz =K p S x (9) 


Sf = T/ K P 
S?j=S/K P 

where T and S are the nominal tensile and shear strengths of 
the interface. If the relative opening displacement S- is not 
zero, the mode mixity can be expressed by 

f>=?f (") 

The mixed-mode damage initiation displacement is 
obtained by substituting Eqs. 7-1 1 into 5, which gives: 


<5° = <?,° 8 


1 + /U 


{S'fHPS'l) 


(K2 


( 12 ) 


The criteria used to predict delamination propagation 
under mixed-mode loading conditions are generally 
established is terms of the energy release rates and fracture 
toughness. The most widely used criteria to predict the 
interaction of the energy release rates in mixed mode is the 
power law given by the expression: 


(gA- 

y G Ic, 


+ 


( r ^ 
G u 

V G Uc j 


(13) 


The exponent a in the power law is usually selected to 
be either 1 or 2, in which case the criterion is a two- 
parameter interaction law with parameters G IC and Guo 

Reeder 22 performed mixed-mode bending (MMB) tests 
to measure the mixed-mode I and II interlaminar fracture 
toughness of composites, providing valuable experimental 
data to assess criteria that have been proposed to predict 
delamination growth. The linear criterion obtained with 
(xr i in Eq. 13 was found to be suitable for thermoplastic 
PEEK matrix composites, its results being comparable to 
more complex interaction functions. However, both the 
linear and quadratic interaction criteria failed to capture 
adequately the dependence of the mixed-mode fracture 
toughness on the mode ratio in epoxy composites. 

In order to accurately account for the variation of 
fracture toughness as a function of mode ratio, a recently 
proposed three -parameter criterion (B-K criterion, Ref. 23) 
is implemented here. The B-K criterion is established in 
terms of the single mode fracture toughnesses G IC and G IIC 
and a parameter Tj : 


g t = G Ic + ( G IIc 




n 


The single-mode failure initiation displacements are then 


(14) 



Reeder 24 applied the B-K criterion to mixed mode test 
results of IM7/977-2 and obtained a best fit using -rj=1.45, 
which is shown in Figure 5. The linear (cc=l) and quadratic 
(oc=2) criteria, which are also shown in Figure 5, do not 
correlate well with the experimental results for this 
material. 



Figure 5. G c versus Gn/G T mode ratio for IM7/99-2. 

For the bilinear stress-displacement softening law 
assumed here, the critical energy release rates in Mode I 
and Mode II are simply the areas under the triangle 0-2-4 in 
Figure 4. 
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J lie 


(15) 


where S F and 8% are the ultimate opening and tangential 
displacements, respectively. 

The Mode I and Mode II energies released at failure 
are computed from: 


Gj 



dSj 



(16) 


where Sj and S 1 ,, are the Mode I and Mode II relative 
displacements at failure under mixed-mode loading. The 
ultimate opening displacement for any mixed-mode 
criterion can be calculated by substituting Eq. 16 into the 
appropriate interaction law. Using the definition of the 
mixed-mode relative displacement S„, in Eq. 8 and the mode 
mixity ratio /i given by Eq. 1 1 , one can solve for the 


ultimate relative displacement S F . For the power law given 
by Eq. 13, the ultimate relative displacement is found to be 


S F 


2(1 + /? 2 ) 
K P S 0 



\ a 


\G, C ) 



{G 1IC J 


l 

a 


(17) 


For the B-K criterion, the ultimate relative displacement is 


8 f =- 


K p S a 


G IC HG 1IC -G IC ) 


p 2 ' 

n~ 

{l+p 2 } 



(18) 


In summary, the mixed mode softening law presented 
above is a single -variable response similar to the bilinear 
single-mode law illustrated in Figure 4. Only one state 
variable, the relative displacement variable 8 m , is used to 
track the damage at the interface. By recording the highest 
value attained by 8m, the unloading response is as shown in 
Figure 4. The displacements for initiation b'°(/j) and 

ultimate failure <J^(/?)are functions of the mode mixity 
P. The relative displacement 8°( P) is computed with Eq. 
12, and s F {p) is computed with either Eq. 17 or Eq. 18. 

The required material parameters are the penalty stiffness 
K p , the interlaminar strengths T and S, the material 
toughnesses G lc and G Uc and either Tj or a. 

The ultimate relative displacements obtained from Eqs. 
17 and 18 can be contrasted to the displacement obtained 
from the mixed mode criterion developed in previous 
work 9 . The criterion in Ref. 9 is based on a quadratic 
interaction between approximations to Eqs. 16, while Eqs. 
1 7 and 1 8 are exact. 

A mixed-mode softening law can be illustrated in a 
single 3D map by representing Mode I on the 2-3 plane, 
and Mode II in the 1-3 plane, as shown in Figure 6. The 
triangle O-T-Sp is the bilinear material response in pure 
Mode I and 0-S-S p is the bilinear material response in 
pure Mode II. It can be observed that the tensile strength T 
is lower than the shear strength S, and the ultimate 
displacement in shear is lower than in tension. In this three- 
dimensional map, any point on the 0-1-2 plane represents a 
mixed-mode relative displacement. 

Under mixed mode, damage initiates at S°, and 
complete fracture is reached at S F . Consequently, the 
tractions for Mode I and Mode II under mixed mode 
loading follow the reduced curves 0-T M -8f M and 
0-S M -S FM , respectively. The areas under these two 
curves represent the fracture energies under mixed mode 
represented by Eqs. 16. 




Figure 6. Combined plot of single-mode bilinear material 
responses. 

The map of all softening responses under mixed mode 
is illustrated in Figure 7. The curve FI represents the 
tractions resulting from the displacements at the onset of 
damage given by Eq. 12, while the curve labeled G 
represents the ultimate relative displacements calculated 
with Eq. 17 or Eq. 18. The triangle 0-A-B is the bilinear 
softening law for a mixed-mode relative displacement of 
S m . The triangle 0-A-B is identical to the triangle 0-2-4 in 
Figure 4. For reference, the triangle 0-C-D in Figure 7 is 
the Mode I bilinear softening response. It can also be 
observed that the effect of compression on the material 
response is neglected. 


Stress 



Figure 7. Map of strain softening response for mixed 
mode delamination. 


three-dimensional (8-node) elements. The relative 
displacements between the top and the bottom faces of the 
element in a local coordinate frame x-y-z are 



where B is the matrix relating the element’s degrees of 
freedom U to the relative displacements between the top 
and bottom interfaces. 

The three-dimensional form of Eq. 3 is 


o = (l-D)C8 or 



( 20 ) 


where I is the identity matrix, C is the undamaged 
constitutive matrix 


’K P 0 0 " 

C = 0 K p 0 (21) 

0 0 K P 

and D is a diagonal matrix representing the damage 
accumulated at the interface: 


~d 0 O' 

D = 0 d 0 (22) 

0 0 d 

The term d on the diagonal is the damage parameter, which 

is a nonlinear function of , the highest mixed-mode 
relative displacement experienced by the material 


S F B max -S° 



Using the maximum value of the relative displacement 
rather than the current value prevents healing of the 
interface. is the only state variable that needs to be 
stored in the database to track the accumulation of damage. 

The minimization of the potential energy subjected to 
the kinematic constraints of Eq. 19 leads to the usual 
integral over the area of the element, which gives the 
following element stiffness 11,21 : 


Element Formulation 

The element stiffness matrix is based on the standard 
isoparametric linear Lagrangian interpolation functions for 


K dta = _[B r ((l-D)c)Bd4 (24) 





The integration is performed numerically using a 
Newton-Cotes integration, which has been shown to 
perform better than Gaussian integration in problems 
involving strain softening 15, lf ’. The integration points of a 
zero-thickness decohesion element coincide with the four 
comers of the element. Since the material softening 
response is evaluated at each integration point, the element 
can soften one comer at a time, giving it the potential to 
model non self-similar delamination growth. 

Nonlinear Solution 

The nonlinear solution of the problems presented here 
was performed using standard ABAQUS procedures. 
However, the softening nature of the interface element 
constitutive equation causes convergence difficulties in the 
solution of the analysis. Schellekens 16 recently suggested 
that in problems where failure is highly localized the 
displacement norm in Riks method should be determined 
considering only the dominant degrees of freedom. May 25 
describes a new automated solution procedure for structures 
with strain-softening materials that is based on a constraint 
equation that uses only the displacement parameters 
associated with the localized failure zone in such structures. 
Unfortunately, a local arc-length procedure was not 
available for the analyses presented here. 

Results and Discussion 

Three test problems were selected to validate the 
decohesion elements. The first problem consists of the 
double cantilever beam (DCB) test used to determine Mode 
I toughness. The second problem modeled consists of the 
end-notched flexure test (ENF) used for Mode II toughness. 
The third test is the mixed bending mode test (MMB). All 
three of these problems have analytical solutions that were 
developed by Mi and Crisfield 19 . These closed form 
solutions provide an approximate framework against which 
to assess the FE models. 

DCB Test for Mode I 

The ASTM standard specimen used to determine the 
interlaminar fracture toughness in Mode I ( (;/<•) is the 
double -cantilever beam (DCB) specimen. This specimen is 
made of a unidirectional fiber-reinforced laminate 
containing a thin insert at the mid-plane near the loaded 
end. A 15-cm.-long specimen, 2 cm.-wide, and composed 
of two 1.98-mm-thick plies of unidirectional material was 
tested by Morais 26 . The initial crack length is 5.5 cm. The 
properties of the graphite material are shown in Table 1, 
and the properties of the interface are shown in Table 2. 


Table 1 . Properties of the Graphite/Epoxy material 26 . 


Ell 

E 22 -E 33 V]2=Vl3 V 23 

C,2=G,3 

G23 

150. GPa 

11. GPa 0.25 0.45 

6.0 GPa 

3.7 GPa 

Table 2 . 

Properties of the DCB specimen interface 26 . 

c, c 

Gnc T 

S 


0.268 N/mm 1.45 N/mm 30. MPa 

40. MPa 

10 6 N/mm 3 


Using Eqs. 9 and the properties of the interface shown 
in Table 2, the relative Mode I and Mode II displacements 
tor damage onset are <5 7 °=30.xl0‘ 6 mm. and 8^ =40.xl0" 6 
mm., respectively. The corresponding ultimate relative 
displacements calculated from Eqs. 15 are Sf = 17.9xl0" 3 
mm. and Sfi =72.5x1 0‘ 3 mm. 

The ABAQUS finite element model, which is shown 
deformed in Figure 8, consists of two layers of C3D8I 
incompatible-mode 8-noded elements. C3D8I elements are 
superior in bending to other low-order continuum elements. 
The anticlastic effects were neglected and only one element 
was used across the width. One hundred and twenty 
elements were used along the span of the model shown in 
Figure 8. 

A plot of reaction force as a function of the applied end 
displacement d is shown in Figure 9. The beam solution 
was developed by Mi and Crisfield 19 for isotropic adherend 
materials and using plane stress assumptions. Note that the 
beam solution is somewhat stiffer than the test and FEM 
results which is probably due to the assumption of isotropy 
in the analytical solution. After the initiation of 
delamination, fiber bridging in the test specimen causes a 
small drift in the response compared to the FEM and 
analytical solutions. 



Figure 8. Model of DCB test specimen. 

Numerical studies with different element sizes indicate 
that the accuracy of the prediction can be significantly 
lower if the size of the elements used in the softening zone 



Applied load 


is greater than some maximum value. The maximum 
predicted load sustained by the DCB specimen calculated 
using several mesh densities is shown in Figure 10. The 
results indicate that poor results are obtained for this 
problem when the element size is greater than 1.25 mm. 
This mesh size is consistent with the results of Goncalves 10 , 
who used 1 -mm. -long 18 -node quadratic elements for the 
analysis of a DCB specimen. 



Figure 11. ENF test specimen. 



Displacement, d, mm. 


Figure 9. Load-deflection response of DCB test. 


The load-deflection responses for the finite element 
model and the analytical prediction are shown in Figure 12. 
It can be observed that both solutions are in excellent 
agreement. 


Mixed Mode Bending 


The most widely used specimen for mixed-mode 
fracture is the mixed-mode bending (MMB) specimen 
shown in Figure 13, which was proposed by Reeder and 
Crews 27 and later re-designed to minimize geometric 
nonlinearities 28 . The mam advantages of the MMB test 
method are the possibility of using virtually the same 
specimen configuration as for Mode I tests and varying the 
mixed mode ratio from pure Mode I to pure Mode II. 



Element Length, mm. Displacement, mm. 


Figure 10. Debond load as a function of element size. Figure 12. Analytical and FEM load-deflection curves for 

ENF test. 


ENF Test for Mode II 


Even though the end-notched flexure test specimen 
shown in Figure 11 exhibits unstable crack propagation for 
short crack lengths, its simplicity makes it a common test to 
measure Mode II fracture toughness. The length of the 
specimen modeled here is 10 cm., its width is 1 cm., and 
the initial debond length is 3 cm. Aluminum adherends 
were used rather than composite to achieve a closer 
approximation to the analytical solutions calculated by 
Mi 19 . The thickness of the adherends is 1.5 mm. The 
properties of the interface are the same as for the DCB 
model. 



Figure 13. MMB test specimen. 

The specimen analyzed here has a length of 10 cm., a 
width of 1 cm., and an initial debond length of 3 cm. The 
thickness of the aluminum adherends is 1.5 mm. The 





properties of the interface are the same for the DCB and 
ENF models. The length of the MMB lever e was chosen 
as 43.72mm, which corresponds to a ratio of 1 for G/G u 
and to a ratio of 2.14 between the load at the mid-span of 
the beam and the opening load. The MMB load fixture is 
simulated by applying an opening load of 100 N. at the 
edge of the debond, and an opposite load of 214 N. at the 
mid-span of the beam. 

The model is composed of two layers of 100 C3D8I 
solid elements. A deformed plot of the finite element model 
is shown in Figure 14. 



Figure 14. Deformed plot of MMB Model. 

The load deflection curve calculated from FEM and 
analytical constant-G curves derived using the beam 
solutions in Ref. 19 are shown in Figure 15 for linear 
interaction between the energy release rates (oc=l), 
quadratic interaction (oe=2), and for a 3 -parameter fit using 
77 = 7 . 45 . In addition, the results of the mixed mode criterion 
from Ref. 9 are also represented. No test results were 
available for the configuration analyzed. However, the 
FEM results compare well with their corresponding 
analytical curves. The FEM results with 77 = 7.45 are based 
on a fit of mixed-mode test data, so they represent the most 
accurate solution and a significant improvement over the 
criterion developed in Ref. 9. 

Concluding Remarks 

A new criterion for the simulation of progressive 
delamination using decohesion elements was presented. 
Decohesion elements are placed between layers of solid 
elements and they open in response to the loading situation. 
The onset of damage and the growth of delamination can be 
simulated without previous knowledge about the location, 
the size, and the direction of propagation of the 
delaminations. A softening law for mixed-mode 
delamination that can be applied to any interaction law was 
proposed. The criterion uses a single state variable, the 
maximum relative displacement , to track the damage 
at the interface under general loading conditions. For the 
linear and quadratic criteria, the material properties 


required to define the element constitutive equation are the 
interlaminar fracture toughnesses and the corresponding 
strengths. The B-K interaction law requires additionally a 
material parameter 77 that is determined from standard 
mixed-mode tests. 

Three examples were presented that test the accuracy 
of the method. Simulations of the DCB and ENF test 
represent cases of single -mode delamination. The MMB 
test that was simulated has equal Mode I and Mode II 
loading conditions. The examples analyzed indicate that the 
mixed-mode criteria can predict the strength of composite 
structures that exhibit progressive delamination. 



Figure 1 5. Predicted load-deflection plots for MMB test. 
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